﻿#pragma once

#include "common.h"
#include <glm/vec2.hpp>
#include <glm/mat3x3.hpp>
#include <glm/mat4x4.hpp>
#include <glm/geometric.hpp>
#include <glm/common.hpp>

namespace anip
{
	using glm::vec2;
	using glm::mat3;
	typedef vec2 invec2; // pass small structs by value
	typedef vec2& outvec2;
	typedef const mat3& inmat3;
	typedef mat3& outmat3;

	wstring toString(invec2 x);

	template<>
	void binaryWrite(const vec2& val, std::ostream& o);

	template<>
	void binaryRead(vec2& val, std::istream& i);

	namespace math
	{
		constexpr double pi = 3.1415926535897932;

		using glm::length;
		using glm::distance;
		using glm::dot;
		using glm::cross;
		using glm::reflect;
		using glm::refract;
		using glm::abs;
		using glm::round;
		using glm::ceil;
		using glm::floor;
		using glm::round;
		using glm::normalize;

		// Note: left-handed, e.g. = (-vy, vx) 
		vec2 normal(invec2 v);

		vec2 rotate(invec2 v, float rad);

		template<typename T>
		T random(T low = 0.0, T high = 1.0)
		{
			return (rand() / (T)RAND_MAX) * (high - low) + low;
		}

		bool inRectangle(invec2 p, invec2 rp, invec2 rs);

		void rectangleIntersect(
			invec2 p1, invec2 s1,
			invec2 p2, invec2 s2,
			outvec2 outp, outvec2 outs);

		template<typename T>
		inline T cubicHermiteInterpolate(float x0, T y0, float x1, T y1,
			T d0, T d1, float x)
		{
			// Given two points on the curve f(x) and their corresponding
			// derivatives as (x0, y0), (x1, y1), y'0 and y'1, a cubic 
			// Hermite interpolation polynomial H(x) for this fragment of 
			// f(x) is a 3rd degree polynomial that satisfies
			//
			//	  H(x0) = y0, H(x1) = y1, H'(x0) = y'0, H'(x1) = y'1.
			//
			// To facilitate solving the coeffecients, we construct H(x) 
			// using four 3rd degree polynomial functions (known as basis 
			// functions) a0, a1, b0 and b1:
			//
			//    H(x) = y0 a0(x) + y1 a1(x) + y'0 b0(x) + y'1 b1(x) [1]
			//
			// where
			//
			//    a0(x0) = 1, a0(x1) = 0, a0'(x0) = 0, a0'(x1) = 0   [2]
			//    a1(x0) = 0, a1(x1) = 1, a1'(x0) = 0, a1'(x1) = 0   [3]
			//    b0(x0) = 0, b0(x1) = 0, b0'(x0) = 1, b0'(x1) = 0   [4]
			//    b1(x0) = 0, b1(x1) = 0, b1'(x0) = 0, b1'(x1) = 1.  [5]
			//
			// Solving [2]-[5] yields
			//
			//                   x - x0    x - x1  2
			//    a0(x) = (1 + 2 -------) (-------)
			//                   x1 - x0   x0 - x1
			//
			//                   x - x1    x - x0  2
			//    a1(x) = (1 + 2 -------) (-------)
			//                   x0 - x1   x1 - x0
			//
			//                      x - x1  2
			//    b0(x) = (x - x0) (-------)
			//                      x0 - x1
			//
			//                      x - x0  2
			//    b1(x) = (x - x1) (-------)
			//                      x1 - x0
			//
			// It can be proved that the H(x) constructed this way is the
			// only Hermite interpolation polynomial under this given
			// condition.
			//    In the following implementation, we assumed x0 == 0 and
			// x1 == 1 for optimization.

			float dist = x1 - x0;
			d0 = d0 / dist; d1 = d1 / dist;
			x = (x - x0) / dist;

			float p1 = powf(1 - x, 2);
			float p2 = x * x;
			float a0 = (1 + 2 * x) * p1;
			float a1 = (3 - 2 * x) * p2;
			float b0 = x * p1;
			float b1 = (x - 1) * p2;
			return y0 * a0 + y1 * a1 + d0 * b0 + d1 * b1;
		}

		template<typename T>
		inline T linearInterpolate(T a, T b, float f)
		{
			return a * (1 - f) + b * f;
		}
	}
}